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Abstract 



We consider the problem of learning a high-dimensional multi-task regression model, un- 
der sparsity constraints induced by presence of grouping structures on the input covariates 
and on the output predictors. This problem is primarily motivated by expression quan- 
titative trait locus (eQTL) mapping, of which the goal is to discover genetic variations 
in the genome (inputs) that influence the expression levels of multiple co-expressed genes 
(outputs), either epistatically, or pleiotropically, or both. A structured input-output lasso 
(SIOL) model based on an intricate £i/£2-norm penalty over the regression coefficient matrix 
is employed to enable discovery of complex sparse input /output relationships; and a highly 
efficient new optimization algorithm called hierarchical group thresholding (HiGT) is de- 
veloped to solve the resultant non-differentiable, non-separable, and ultra high-dimensional 
optimization problem. We show on both simulation and on a yeast eQTL dataset that our 
model leads to significantly better recovery of the structured sparse relationships between 
the inputs and the outputs, and our algorithm significantly outperforms other optimization 
techniques under the same model. Additionally, we propose a novel approach for efficiently 
and effectively detecting input interactions by exploiting the prior knowledge available from 
biological experiments. 

Keywords: Group Lasso, Multi-task Lasso, Epistasis, Pleiotrophy, Genome-wide associa- 
tion studies, eQTL mapping, Genetic interaction network 
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1. INTRODUCTION 
In this paper, we consider the problem of learning a functional mapping from a high- 
dimensional input space with structured sparsity to a multivariate output space where re- 
sponses are coupled (therefore making the estimator doubly structured), with an application 
for detecting genomic loci affecting gene expression levels, a problem known as expression 
quantitative trait loci (eQTL) mapping. In particular, we are interested in exploiting the 
structural information on both the input and output space jointly to improve the accuracy for 
identifying a small number of input variables relevant to the outputs among a large number 
of candidates. When the input or output variables are highly correlated among themselves, 
multiple related inputs may synergistically influence the same outputs, and multiple related 
outputs may be synergistically influenced by the same inputs. 

The primary motivation for our work comes from the problem of genome- wide association 



(GWA) mapping of eQTLs in computational genomics [2l|, of which the goal is to detect the 
genetic variations, often single nucleotide polymorphisms (SNPs), across the whole genome 
that perturb the expression levels of genes, given the data of SNP genotypes and microarray 
gene-expression traits of a study cohort. One of the main challenges of this problem is 
that, typically the sample size is very small (e.g. ~ 1000), whereas there are a very large 
number of SNPs (e.g. ~ 500,000) and expression traits (e.g., ~ 10,000). Furthermore, 
there have been numerous evidences that multiple genetic variations may interact with each 



other (a.k.a., epistasis) 



39 



genes (a.k.a., pleiotropy) [41 



ll, and the same genetic variation(s) can influence multiple 
14j . However, prior knowledge of such information implying 



complex association structures are difficult to exploit in standard statistical analysis of GWA 



mapping [48 



321 ] . To enhance the statistical power for mapping of eQTLs, it is desirable 
to incorporate biological knowledge of genome and transcriptome structures into the model 
to guide the search for true eQTLs. In this article, we focus on developing a model which 
can make use of structural information on both input (SNPs) and output (gene expressions) 
sides. In particular, we consider biological knowledge about group associations as structural 
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information. If there exist group behaviors among the covariates in the high-dimensional 
input X, for example, multiple genetically coupled SNPs (e.g., in linkage disequilibrium) 



can jointly affect a single trait 



48], such group information is called an input structure; 



if multiple variables in the high-dimensional output Y are jointly under the influence of a 
similar set of input covariates, for example, a single SNP can affect multiple functionally 



23], such group information is 



coupled traits (e.g., genes in the same pathway or operon) 
called an output structure. The problem of GWA mapping of eQTLs can be formulated 
as a model selection problem under a multitask regression model Y = BX with structured 
sparsity, where the resultant non-zero elements in the regression coefficient matrix B expose 
the identities of the eQTLs and their associated traits. 

Variants of this problem have been widely studied in the recent high-dimensional infer- 
ence and variable selection literature, and various penalty-based or Bayesian approaches for 



learning a shared sparsity pattern amon g ei t 
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rer multiple inputs or multiple outputs in a 
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regression model have been proposed [45 
constraints, different penalty functions have been previous 
norm (group-lasso) penalty for a simple grouping structure 



271 ] . Depending on the type of structural 
y considered, including mixed- 



52 



penalty for a tree structure 



54] . tree-guided group-lasso 



22] , or graph-guided fused lasso for a graph structure [23[ . Most 



previous approaches, however, considered either only the input structural constraints or only 
the output structural constraints, but not both. There have been a few approaches that at- 



tempted to use both structural information, including MGOMP 



multi-task learning" 



and "group sparsity for 



4=71 ] . MGOMP proposed to select the groups of regression coefficients 



471 ] proposed to find the 



from a predefined set of grouped variables in a greedy fashion, and 
groups of inputs that influence the group of outputs. However, both methods may have lim- 
its on the number or the shapes of sparsity patterns that can be induced in B. For example, 
given a large number of input groups \Q\ and output groups \H\ (e.g. \Q\ > 10 5 , \7i\ > 10 3 
for genome data) the scalability of MGOMP can be substantially affected since it needs to 
select the groups of coefficients from all possible combinations of input and output groups. 



For [47] . only disjoint block sparsity patterns are considered, hence it may not capture the 



sparsity patterns where the grouped variables overlap. 

In this paper, we address the problem of exploiting both the input and output structures 
in a high-dimensional linear regression setting practically encountered in eQTL mapping. 
Furthermore, to detect epistatic (i.e., interaction) effects between SNP pairs, we addition- 
ally expand the input space to include pairwise terms ( 's) guided by biological 
information, which necessitates attentions for avoiding excessive input dimension that can 
make the problem computationally prohibitive. Our main contributions can be summarized 
as follows: 

1. We propose a highly general regression model with structured input-output regular- 
izes called "jointly structured input-output lasso" (SIOL) that discovers structured 
associations between SNPs and expression traits (Section [3]). 

2. We develop a simple and highly efficient optimization method called "hierarchical 
group thresholding" (HiGT) for solving the proposed regression problem under com- 
plex sparsity- inducing penalties in a very high-dimensional space (Section H]). 

3. Extending SIOL, we propose "structured polynomial multi-task regression" to effi- 
ciently model non-additive SNP-SNP interactions guided by genetic interaction net- 
works (Section [5]). 

Specifically, given knowledge of the groupings of the inputs (i.e., SNPs) and outputs (i.e., 
traits) in a high-dimensional multi-task regression setting, we employ an L1/L2 norm over 
such structures to impose a group-level sparsity-inducing penalty simultaneously over both 
the columns and the rows of the regression coefficient matrix B (In our setting, a row cor- 
responds to coefficients regressing all SNP (or SNP pair) inputs to a particular trait output, 
thus reflecting possible epistatic effects; and a column corresponds to coefficients regressing 
a particular SNP (or SNP pair) input to all trait outputs in question, thus reflecting possible 
pleotropic effects). Given reliable input and output structures, rich structured sparsity can 
increase statistical power significantly since it makes it possible to borrow information not 
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only within different output or input variables, but also across output and input variables. 

The sparsity-inducing penalties on both the inputs and outputs in SIOL introduce a non- 
different iable and non-separable objective in an extremely high-dimensional optimization 



space, which prevents standard optimization methods such as the interior point 351 ]. the 



coordinate-descent [15J], or even the recently invented union of supports [19] algorithms 
to be directly applied. We propose a simple and efficient algorithm called "hierarchical 
group-thresholding" to optimize our regression model with complex structured regularizers. 
Our method is an iterative optimization algorithm, designed to handle complex structured 
regularizers for very large scale problems. It starts with a non-zero B (e.g. initialized by ridge 



regression 18] ) , and progressively discards irrelevant groups of covariates using thresholding 
operations. In each iteration, we also update the coefficients of the remaining covariates. To 
speed up our method, we employ a directed acyclic graph (DAG) where nodes represent the 
zero patterns encoded by our input-output structured regularizers at different granularity, 
and edges indicate the inclusion relations among them. Guided by the DAG, we could 
efficiently discard irrelevant covariates. 

As our third contribution, we consider non-additive pairwise interaction effects between 
the input variables, in a way that avoids a quadratic blow-up of the input dimension. In 
eQTL mapping studies, it is not uncommon that the effect of one SNP on the expression-level 
of a gene is dependent on the genotype of another SNP, and this phenomenon is known as 
epistasis. To capture pairwise epistatic effects of SNPs on the trait variation, we additionally 
consider non-additive interactions between the input covariates. However, in a typical eQTL 
mapping, as the input lies in a very high- dimensional space, it is computationally and statis- 
tically infeasible to consider all possible input pairs. For example, for J inputs (e.g. 500, 000 
for a typical genome data set), we have 0(J 2 ) candidate input pairs, and learning with all 
of them will require a significantly large sample size. Many of the previous approaches for 
learning the epistatic interactions relied on pruning candidate pairs based on the observed 



data 



38 



or constructing candidate pairs from individual SNPs that were selected 



Dased on 



marginal effects in the previous learning phase without modeling interactions 49 



13|. A 



main disadvantage of the later approach is that it will miss pairwise interactions when they 
have no or little individual effects on outputs. Instead of choosing candidate SNP pairs based 



ldj constructed from 



on only marginal effects, we propose to use genetic interaction network 
large-scale biological experiments to consider biologically plausible candidate pairs. 

The rest of this paper is organized as follows. In Section 2, we discuss previous works 
on learning a sparse regression model with prior knowledge on either output or input struc- 
ture. In Section 3, we introduce our proposed model "jointly structured input-output lasso" 
(SIOL). To solve our regression problem, we present an efficient optimization method called 
"hierarchical group-thresholding" (HiGT) in Section 4. We further extend our model to 
consider pairwise interactions among input variables and propose "structured polynomial 
multi-task regression" in Section 5. We demonstrate the accuracy of recovered structured 
sparsity and the speed of our optimization method in Section 6 via simulation study, and 
present eQTLs having marginal and interaction effects in yeast that we identified in Section 
7. A discussion is followed in Section 8. 



2. BACKGROUND: LINEAR REGRESSION WITH STRUCTURED SPARSITY 
In this section, we lay out the notation and then review existing sparse regression methods 
that recover a structured sparsity pattern in the estimated regression coefficients given prior 
knowledge on input or output structure. 

2.1 Notation for matrix operations 

Given a matrix B £ M. KxJ , we denote the fc-th row by (3^, the j-th column by f3 J , and the 
(k,j) element by (3 3 k . \\-\\f denotes the matrix Frobenius norm, ||-||i denotes an Li norm 
(entry- wise matrix L\ norm for a matrix argument), and ||-||2 represents an L 2 norm. Given 
the set of column groups Q = {gi, . . . , g\g\} defined as a subset of the power set of {1, ... , J}, 
f3f represents the row vector with elements {/3 3 k : j ; £ g, g £ Q}, which is a subvector of f3k 
due to group g. Similarly, for the set of row groups % = {hi, . . . , h^} over M rows of matrix 
B, we denote by (3 3 h the column subvector with elements {f5 k : k £ h, h £ H}. We also define 
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the submatrix of as a |h| x |g| matrix with elements {(3 J k : k G h, j G g, h G H, g G Q}. 
2.2 Sparse estimation of linear regression 

Let X G IR Jx7V be the input data for J inputs and N individuals, and Y G IR Kx7V be the out- 
put data for K outputs. We model the functional mapping from the common J-dimensional 
input space to the K- dimensional output space, using a linear model parametrized by un- 
known regression coefficients B G ]R^ xJ as follows: 

Y = BX + E, 

where E G M K xN is a matrix of noise terms whose elements are assumed to be identically and 
independently distributed as Gaussian with zero mean and the identity covariance matrix. 
Throughout the paper, we assume that x*'s and yj^'s are standardized such that all rows of X 
and Y have zero mean and a constant variance, and consider a model without an intercept. 
In eQTL analysis, inputs are genotypes for J loci encoded as 0, 1, or 2 in terms of the 
number of minor alleles at a given locus, and output data are given as expression levels of 
genes measured in a microarray experiment. Then, the regression coefficients represent the 
strengths of associations between genetic variations and gene expression levels. 

Our proposed method for estimating the coefficients B is based on a group-structured 
multi-task re gres sion approach that extends existing re gula rized regression approaches in- 



cluding lasso |44|, group lasso [52| and multi-task lasso [36|, which we briefly review below 
in the context of our eQTL mapping problem. When J » N and only a small number of 
inputs are expected to influence outputs, lasso has been widely used and shown effective in 
selecting the input variables relevant to outputs and setting the elements of B for irrelevant 



inputs to zero 



53]. Lasso obtains a sparse estimate of regression coefficients by optimizing 



the least squared error criterion with an L\ penalty over B as follows: 

imn^||Y-BX|||. + A||B|| 1) (1) 

where A is the tuning parameter that determines the amount of penalization. The optimal 
value of A can be determined by cross validation or via an information-theoretic test based on 



BIC. As in eQTL analysis it is often believed that the expression level of each gene is affected 
by a relatively small number of genetic variations in the whole genome, lasso provides an 
effective tool for identifying eQTLs from a large number of genetic variations. Lasso has 
been previously applied to eQTL analysis |5| and more general genetic association mapping 
problems 50 1. 

While lasso considers the input variables independently to select relevant inputs with 
non-zero regression coefficients, we may have prior knowledge on how related input variables 
are grouped together and want to perform variable selection at the group level rather than at 
the level of individual inputs. Grouped variable selection approach can combine the statis- 
tical strengths across multiple related input variables to achieve higher power for detecting 
relevant inputs in the case of low signal-to-noise ratio. Assuming the grouping structure over 
inputs are available as Q = {gi, . . . , g\g\}, which is a subset of the power set of {1, . . . , J}, 
group lasso uses Li/L 2 penalization to enforce that all of the members in each group of input 
variables are jointly relevant or irrelevant to each output. Group lasso obtains an estimate 



of B by solving the following optimization problem: 

K 

IIY-BXIIS. + A' 

b 2 



mini||Y-BX||| + A^^||/3f|| 2 , (2) 

k=i gee 

where A is the tuning parameter. The second term in the above equation represents an 
Li/L 2 penalty over each row (3u of B for the k-th. output given Q, defined by ||/3fc||L 1 /L 2 = 
J2 s egW@k lb- The L 2 part of the penalty plays the role of enforcing a joint selection of inputs 
within each group, whereas the L\ part of the penalty is applied across different groups to 
encourage a group-level sparsity. Group lasso can be applied to an eQTL mapping problem 
given biologically meaningful groups of genetic variations that are functionally related. For 
example, rather than individual genetic variations acting independently to affect (or not 
affect) gene expressions, the variations are often related through pathways that consist of 
multiple genes participating in a common function. Thus, genetic variations can be grouped 
according to pathways that contain genes carrying those genetic variations. Then, given this 
grouping, group lasso can be used to select groups of genetic variations in the same pathways 
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as factors influencing gene expression levels 40] 



Instead of having groups over inputs with outputs being independent as in group lasso, 
the idea of using L\jL 2 penalty for grouped variable selection has also been applied to take 
advantage of the relatedness among outputs in multiple output regression. In multi-task 
regression for union support recovery 3f| , one assumes that all the outputs share a common 
support of relevant input variables and try to recover shared sparsity patterns across multiple 
outputs by solving the following optimization problem: 

1 J 
mm-||Y-BX||^ + A^^||^ h || 2 , (3) 

j=\ hew 

where A can be determined by cross-validation. In eQTL mapping, as gene expression levels 

are often correlated for the genes that participate in a common function, it is reasonable 

to assume that those coexpressed genes may be influenced by common genetic variations. 

If gene module information is available, one can use the above model to detect genetic 

variations influencing the expressions of a subset of genes within each gene module. This 

strategy corresponds to a variation of the standard group lasso, where group is defined over 

outputs rather than inputs. 

Extending the idea of lasso and group lasso, we may have group and individual level 

sparsity simultaneously using combined L\ and Li/L 2 penalty. In group lasso, if a group 

of coefficients is not jointly set to zero, all the members in the group should have non-zero 

values. However, sometimes it is desirable to set some members of the group to zero if they 

are irrelevant to outputs. Sparse group lasso 

groups of coefficients include both relevant and irrelevant ones. Using convex combination 

of Li and Li/L 2 norms, it solves the following convex optimization problem: 

i K 
mm - ||Y - BX||J. + X 1 \\B\\ 1 + A 2 £ £ WW* (4) 

k=i gee 

where Ai and X 2 determine the individual and group level sparsity, respectively. The L\jL 2 
penalty shrinks groups of coefficients to zero, and at the same time, L\ penalty sets irrelevant 
coefficients to zero individually within each group. 

10 



161 ] is proposed to address the cases where 



Our proposed model is motivated by group lasso, multi-task lasso and sparse group lasso, 
each of which can exploit pre-defined groupingness of input or output variables to achieve 
better statistical power. In the next section, we will extend the existing models in such a 
way that we can use the groups in both input and output spaces simultaneously. Adopting 
the idea of sparse group lasso, we will also support variable selection at individual levels. 

3. JOINTLY STRUCTURED INPUT-OUTPUT LASSO 
In this section, we propose SIOL that incorporates structural constraints on both the inputs 
and outputs. The model combines the mixed-norm regularizers for the groups of inputs and 
outputs, which leads to the following optimization problem: 

min-||Y-BX||^ + Ai||B||i, (5a) 

+ ^EEn» ( 5b ) 

k=i gee 

+ As££ll#IU (5c) 

j=l h£H 

where Eq. (l5b|) incorporates the groups of inputs Q = {gi, . . . , g\g\}, Eq. ( |5c|) incorporates 
the groups of the outputs H = {hi, . . . , h^}, and Eq. (15a]) allows us to select individual 
coefficients. Note that it is possible that there are overlaps between /3f and f3f , between 
f3 3 h and j3 ] h ,, and between f3f and f3 3 h , where g ^ g', h ^ h' and g, g' G Q, h, h' G H. The 
overlaps make it challenging to optimize Eq. ((5]), and this issue will be addressed by our 
optimization method in Section |H 

Let us characterize the structural constraints imposed by the penalties in our model. In 
our analysis, we investigate a block of coefficients involved in one output group h and one 
input group g, i.e, B^. We start with Karush-Kuhn- Tucker (KKT) condition for Eq. ([5]): 

(y k - A-X)( Xj ) T = Ai4 + A 2 c£ + \ 3 d{, (6) 

where s 3 k , c? k , and d\ are the subgradient of Eq. fl5a|) . Eq. f lFbj) . and Eq. fl5c|) with respect to 
01, respectively. For simple notation, we also define r J k = — /4 X «- 
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First, we consider the case where all coefficients in become zero simultaneously, i.e., 
B| = 0. Using KKT condition in Eq. fl6]), we can see that = if and only if 

EEI r iN r - M} 2 < EE (M + M) 2 < (a 2V ^T+ a 3V ^) 2 • (r) 

feeh jeg fceh jeg 

This condition is due to Cauchy-Schwarz inequality, ^ jGg (Cfc) 2 < 1, and Xlfceh(^fc) 2 — 1- 
Here if Ai, A2 and A3 are large, B^ is likely to be zero jointly. This structural sparsity is 
useful to filter out a large number of irrelevant covariates since it considers both the group 
of correlated inputs g and the group of correlated outputs h simultaneously. 

Our model also inherits grouping effects for only input (or output) groups. For the 
analysis of such grouping effects, we fix the groups of zero coefficients that overlap with, 
say, an input group (3f. Formally speaking, let us define £ = {j : (f3 ] h , = 0,j 6 g, h' G 
%) v (/3f' = 0, j G g' A g)}, and fix f3{s for all j G £. Using the KKT condition in Eq. (ED, 

ft = if 

E ~ x ^ ^ E + A 3 rf i) 2 < 4 (8) 

Here, we know that d{ = for j G g - £ ($( = and ^ 0) and A 2 £ ieg (/?£) 2 = 
A2y^„>-„ c(0t) 2 , and hence ^ gg _£ (A2<4 + A 3 c^) < A|. This technique was previously 
to handle overlapping group lasso. One can see that if the size of £ is 
large, (3f tends to be zero together since it reduces the left-hand side of Eq. (jSJ). This behavior 
explains the correlation effects between input and output group structures. When a group of 
coefficients (J3f, f3 J h ) corresponding to an input group or an output group become zero, they 
affect other groups of coefficients that overlap with them; and the overlapped coefficients are 
more likely to be zero. These correlation effects between overlapping groups are desirable for 
inducing appropriate structured sparsity as it allows us to share information across different 
inputs and different outputs simultaneously. We skip the analysis of the grouping effects 
for output groups as the argument is the same except that the input and output group are 
reversed. 



introduced in 51 
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Finally, we also have individual sparsity due to L\ penalty in Eq. (15 ap . In this case, let 
us assume that (3f 7^ and {3^ ^ since if the group of coefficients is zero, we automatically 
have j3 3 k = 0. Using the KKT condition, (3^ = if and only if 

|ri(x,) T | < Ax. (9) 

It is equivalent to the condition of lasso that sets a regression coefficient to zero. Note that 
if A2 = A3 = 0, we have sparsity only at the individual levels, and our model is the same as 
lasso. When a group of coefficients contains both relevant and irrelevant ones, we can set 
the irrelevant coefficients to zero using Eq. fl9]). 

We briefly mention the three tuning parameters (Ai, A2, A3) which can be determined by 
cross validation. It is often computationally expensive to search for optimal parameters in 
3-dimensional grid. In practice, instead, we use the following tuning parameters: A2 = A 2 Ag 
and A3 = (1 — A 2 )A 3 . Here A2 configures the mixing proportion of input and output group 
structures, and A3 is the scaling factor that determines the degree of penalization for the 
input and output groups. In this setting, we also have three regularization parameters, 
however, it helps us to reduce the search space of the tuning parameters as we know the 
range of A' 2 (0 < A 2 < 1). 

Let us discuss the statistical and biological benefits of our model. First, our model 
can capture rich structured sparsity in B. The structured sparsity patterns include zero 
(sub)rows, zero (sub)columns and zero blocks of B^. It is impossible to have such rich 
sparsity patterns if we use one part of information on either input or output side. For 



example, group lasso 52[ or multi-task lasso 36[ consider structured sparsity patterns in 
either rows or columns in B. Second, our model is robust to the groups which contain both 
relevant and irrelevant coefficients. If predefined groups of inputs and outputs are unreliable, 
our model may still work since the irrelevant coefficients can be set to zero individually via 
Li penalty even when their groups are not jointly set to zero. Third, the grouping effects 
induced by our model in Eq. ([TJ [8]) show that we can use the correlation effects between 
input and output groups. When we have reliable input and output groups, the advantage 
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from the structural information will be further enhanced by the correlation effects in addition 
to the sum of the benefits of both input and output groups. 

When applied to GWA mapping of eQTLs, our model offers a number of desirable prop- 
erties. It is likely that our model can detect association SNPs with low signal-to-noise ratio 
by taking advantage of rich structural information. In GWA studies, one of the main chal- 
lenges is to detect SNPs having weak signals with limited sample size. In complex diseases 
such as cancer and diabetes, biologists believe that multiple SNPs are jointly responsible for 



diseases but not necessarily with strong marginal effects [30[. However, such causal SNPs are 
hard to detect mainly due to insufficient number of samples. Our model can deal with this 
challenge by taking advantage of both input and output group structures. First, by grouping 
inputs (or SNPs), we can increase the signal-to-noise ratio. Suppose each SNP has small 
signal marginally, if a group of coefficients is relevant, their joint strength will be increased, 
and it is unlikely that they are jointly set to zero. On the other hand, if a group of coeffi- 
cients is irrelevant, their joint strength will still be small, and it is likely that they are set to 
zero. Second, taking advantage of the output groups, we can share information across the 
correlated outputs, and it decreases the sample size required for successful support recovery 



341 ] . Overall, to detect causal SNPs having small effects, our model increases signal-to- noise 
ratio by grouping the SNPs, and simultaneously decreases the required number of samples 
by grouping phenotypic traits. 

Unfortunately, the optimization problem resultant from Eq. (JSJ) is non-trivial. One 
may find out that each (3 J k appears in all the three penalties of Eq. fl5al - l5cl). Thus, our 
structured regularizer is non-separable, which makes simple coordinate descent optimization 
inapplicable. The overlaps bet ween/ within input and output groups add another difficulty. 
Furthermore, we must induce appropriate sparsity patterns (i.e., exact zeros) in addition to 
the minimization of Eq. (JHJ), therefore approximate methods based on merely relaxing the 
shrinkage functions are not appropriate. In the following section, we propose "hierarchical 
group thresholding" method (HiGT) that efficiently solves our optimization problem with 
hierarchically organized thresholding operations. 

14 



(a) 




(b) 

Figure 1: Sparsity patterns of and a DAG constructed with the sparsity patterns. The 
shaded area shows zero entries, (a) All possible zero patterns of that can be induced by 
Eq. f)5a|5b|5cp when g = {1,2} and h = {1,2}. (b) An example of a DAG that contains 
the zero patterns of B^. The root node contains zero pattern for B^ = 0, and the internal 
nodes represent the zero patterns for f3 J h = (one column is zero) or f3f = (one row is 
zero). The leaf nodes denote fi 3 k = 0. In the DAG, the zero pattern of children nodes should 
be a subset of their parent nodes' zero patterns. 



4. OPTIMIZATION METHOD 
In this section, we propose our method to optimize Eq. ([5]). We start with a non-zero B 
initialized by other methods (e.g. ridge regression), and always reduce the set of non-zero 
/3^s using thresholding operations as our procedure proceeds. Our framework is an iterative 
procedure consisting of two steps. First, we set the groups (or individual) of regression 
coefficients to zero by checking optimality conditions (called thresholding) as we walk through 
a predefined directed acyclic graph (DAG). When we walk though the nodes in the DAG, 
some (3 3 k s might not achieve zero. Second, we update only these non-zero f3 3 k s using any 
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available optimization techniques. 

Let us first characterize the zero patterns induced by Eq. (I3al - 13c]) . We separately 
consider a block of B which consists of one input group (g G Q) and one output group 
(h G "H). Our observation tells us that there are grouping effects (to be zero simultaneously) 
for each g and h: /3f = and (3^ = 0. We also have = when /3f = 0, VA; G h 
or {3 J h = 0, Vj G g. Each covariate can also be zero, i.e, /3 3 k = due to the l\ penalty 
in Eq. f)5ap . Figure [IJa) shows all the possible zero patterns of B^ induced by Eq. (IBlil 
- l5cl). Given these sparsity patterns, to induce structured sparsity, one might be able to 
check whether or not these zero patterns satisfy optimality conditions and discard irrelevant 
covariates accordingly. However, this approach may be inefficient as it needs to examine 
the large number of zero patterns. Instead, to efficiently check the zero patterns, we will 
construct a DAG, and exploit the inclusion relationships between the zero patterns. The 
main idea is that we want to be able to check all zero patterns by traversing the DAG while 
avoiding unnecessary optimality checks. 

In Figure [Hb), we show an example of the DAG for B^ when g = {1, 2} and h = {1, 2}. 
We denote the set of all possible zero patterns of B by Z = {Zi, . . . , Z\z\}- For example, 
Z\ can be a zero pattern for B^ = (the root node in Figure G^b)). Let us denote B(Z t ) 
by the coefficients of B corresponding to ZfS zero pattern. Then we define the DAG as 
follows: A node is represented by Z G Z, and there exists a directed edge from Z\ G Z to 
Z 2 G Z if and only if Z\ D Z 2 and $Z G Z : Z\ D Z D Z 2 . For example, in Figure H^b), the 
zero patterns of the nodes in the second level include the zero patterns of their children. In 
general, when we have multiple input and output groups, we can generate a DAG for each 
B| separately and then connect all the DAGs to the root node for B = 0. This graph is 
originated from Hasse diagram 6] and it was previously utilized for finding a minimal set of 



groups for inducing structured sparsity 20 1. 

We can readily observe that our procedure has the following properties: 
• Walking through the DAG, we can check all possible zero patterns explicitly without 
resorting to heuristics or approximations. 
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• If B(Z) = 0, Z £ Z, we know that all the descendants of Z are also zero due to the 
inclusion relations of the DAG. Hence, we can "skip" to check the optimality conditions 
that the descendants of Z are zero. 

Considering these properties, we develop our optimization framework for the following 
reasons. First, we can achieve accurate zero patterns in B since we check all possible sparsity 
patterns when walking through the DAG. Second, if B is sparse, our framework is very 
efficient since we can skip the optimality checks for many zero patterns in Z. Mostly we will 
check only nodes located at the high levels of the DAG. Third, our framework is simple to 
implement. All we need is to check whether each node in the DAG attains zero and update 
non-zero j3 3 k s only when necessary. 

Specifically, our hierarchical group-thresholding employs the following procedure: 

1. Initialize a non-zero B using any available methods (e.g. ridge regression). 

2. Construct a DAG that contains all zero patterns of B that can be induced by the 
penalty in Eq. ( T5al l5bl l5cl) . 

3. Use depth-first-search (DFS) to traverse the DAG, and check the optimality conditions 
to see if the zero patterns at each node Z achieve zero. If B(Z) = or Z satisfies the 
optimality condition to be zero, set B(Z) = 0, skip the descendants of Z, and visit the 
next node according to the DFS order. 

4. For those of /3^'s which did not achieve zero in the previous step, update the coefficients 
of the non-zero /3^'s using any available optimization algorithms. 

5. Iterate step [3] and H] until Eq. ([5]) converges. 

Bellow we briefly present the derivations of the three ingredients of our optimization 
framework that include 1) the construction of a DAG, 2) the optimality condition of each 
Z G Z in the DAG and 3) the rule for updating non-zero regression coefficients. Our 
optimization method is summarized in Algorithm [TJ 
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Algorithm 1 Hierarchical group-thresholding method for Eq. (jSJ) 

B <— coefficients estimated by ridge regression 
Q <— groups of inputs 
"H <— groups of outputs 

D(Z,£) <— DAG including all zero patterns 
{Z (1) ,Z (2) , . <- DFS order of 2 in D 

repeat 
i <- 1 

while t < \Z\ do 

if contains = then 

c <r- Eq. ITOt 
else if Zm contains /3f = then 

c <- Eq. JTT1 
else if Z/ t \ contains 0^ = then 

c <r- Eq. lT2t 
else if Zha contains fP h = then 

c <- Eq. lT3t 
end if 

if c holds (condition for B(Z( t j) = 0) or B(Zm) = then 
B(Z( t j) = (Set zero to Z(t)' s zero pattern) 

t «— DFS order of t' such that Z( t /j is not a descendant of Z/ t \ 7 t' > t and jit" : t' > t" > t (Skip the descendants 

of z (t)) 
else if c = Eq. 11131 then 

Update f3 3 k using Eq. i l 1 4[ ) (Updating non-zero regression coefficients) 

t <r- t + 1 

else 

t «- t + 1 
end if 
end while 
until convergence 

Construction of the DAG To generate the DAG, first we define the set of nodes Z. For 
each block of B n , we are interested in the four types of zero patterns as follows: 

1. B n is zero: = 0. 

2. One row in B^ is zero: /3| = 0, k G h. 

3. One column in B^ is zero: f3 J h = 0, j G g. 

4. One regression coefficient in B n is zero: f3 3 k = 0, fc G h and j G g. 
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These zero patterns of B are shown in Figure[]Jb) when |g| = |h| = 2. For example, Case 
2 and 3 correspond to the nodes at the second level of the DAG. For all g G Q and h G H, 
we can define nodes Z G Z using the above zero patterns. Then we need to determine the 
edges of the DAG by investigating the relations of the nodes. We can also easily see that 
there exists the relationship among the zero patterns: Case 1 D Case 2, Case 3 D Case 4. 
Given the zero patterns and their relations, we create a directed edge Zi — > Z 2 if and only 
if Z\ D Z 2 and pZ G Z : Z\ Z> Z D Z 2 - In Figure QJb) we show an example of the DAG. 
Finally, we make a dummy root node and generate an edge from the dummy node to the 
root of all DAGs for B| = 0. 

Optimality conditions for structured sparsity patterns Given a block of B^, here 
we show optimality conditions for the four sparsity patterns: (1) B^ = 0, (2) f3f = 0, (3) 
(3 3 h = 0, and (4) J k = 0, (j G g, k G h, g G Q, h G H). In Figure d^b), the root node 
corresponds to the first case, the nodes at the second level correspond to the second and 
third case, and the leaf nodes correspond to the fourth case. Our derivation of the following 
optimality conditions use the fact that all zero coefficients are fixed, as it makes it simple 
to deal with overlapping groups. We denote the column and row indices of zero entries by 
r) = {j : (3{ = 0, Vj G g, V&; G h} and 7 = {k : /3 J k = 0, Vj G g, VA; G h}. 
First, the optimality condition for the first case is as follows: B^ = if 



E E Hte) T - A ^} 2 < + A 3V ^) , (10) 



where 



fceh-7 jeg-r) 



Ai 



if 



' r J ( x ) T 

sign I k \ ) if 



ri(x,) 5 



< 1 
> 1. 



It is derived using KKT condition in Eq. ([6]) and Cauchy-Schwarz inequality. 
The second case of structured sparsity, i.e, /3f = is achieved if 



E H^r - \ lS i} 2 < xi, ai) 
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and the optimality condition for the third case, i.e, (3 3 h = is 



£ {4(x,-) T - A^} 2 < A 
fceh — / 



121 



These conditions can be established using KKT condition in Eq. fl6]) fixing all the zero 
coefficients. Finally, assuming that f3 J h ^ and f3f ^ 0, the fourth case has the optimality 
condition of 



|ri(x,fl 



(13) 



Update rule for nonzero coefficients If all the above optimality conditions do not hold, 
we know that (3 k ^ 0. In this case, the gradient of Eq. ([5]) with respect to (i k exists, and 
we can update {3 3 k using any coordinate descent procedures. With a little bit of algebra, we 
derive the following update rule: f3 k = j3 3 k _ + (3 k + where 



hi 



mm 



max 



E 



A, 



*S 11^ 



o, i+j: 



^ H^fclla ki£ 11^ 



l 

{ratef-Ai} 



(14) 



We close this section by summarizing the desirable properties of our optimization method. 
First, when B is sparse, our optimization procedure is very fast. We take advantage of not 
only the hierarchical structure of the DAG, but also the simple forms of the optimality 
conditions with residuals. If we keep track of the residuals, we can efficiently check the 
optimality conditions for each sparsity pattern. Second, our thresholding operations check all 
possible sparsity patterns, resulting in appropriate structured sparsity in B. It is important 
for eQTL mapping since the coefficients for irrelevant SNPs can be set to exactly zero. Third, 
our optimization method can deal with overlaps between/within the coefficients for input 
groups (/3f's) and output groups (/3h's). Since input or output groups may overlap, and they 
must be considered simultaneously, this property of our method is essential. Finally, unlike 



some previous methods 



52 



44|, we make no use of the assumption that the design matrix X 
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is orthonormal (X T X = I). This dropping of the assumption is desirable for eQTL mapping 
in particular as covariates (SNPs) are highly correlated due to linkage disequilibrium. If 
one uses orthonormalization as a preprocessing step to make X orthonormal, there is no 



guarantee that the same solution for the original problem is attained 16]. 



5. DEALING WITH STRUCTURES INDUCING HIGHER-ORDER EFFECTS 
So far, we have been dealing with input and output structures in the context of multi- 
variate and multi-task linear regression where the influences from the covariates on the 
responses are additive. When higher interactions take place among covariates, which is 



known as epistasis and is prevalent in gen etic associations 



such effects is polynomial regression 



7j, a common approach to model 



31| . where higher-order terms of the covariates are 
included as additional regressors. However, in high-dimensional problems such as the one 
studied in this paper, this strategy is infeasible even for 2nd-order polynomial regression 
because, given say, even a standard genome dataset with ~ 10 5 SNPs, one is left with ~ 10 
regressors which is both computationally and statistically unmanageable. In this section, 
we briefly show how to circumvent this difficulty using structured regularization based on 
prior information of covariate interactions. This strategy is essentially a straightforward 
generalization of the ideas in Section [3] to a polynomial regression setting using a special 
type of structure encoded by a graph. Therefore all the algorithmic solutions developed in 
Section H] for the general optimization problem in Section [3] still apply here. 

Following common practice in GWA literature, here we consider only 2nd-order inter- 
actions between SNP pairs. Instead of including all SNP pairs as regressors, we employ 
a synthetic genetic interaction network 



ldj to define a relatively small candidate set U of 
interacting SNP pairs. A synthetic genetic interaction network is derived from biological 
evidences of pairwise functional interactions between genes, such as double knockout experi- 



ments 



46 



10 



, y| . It contains information about the pairs of genes whose mutations affect 
the phenotype only when the mutations are present on both genes, and this represents a set 
of ground-truth interaction effects. Given such a network, we consider only those pairs of 
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SNPs that are physically located in the genome near the genes that interact in the network 
within a certain distance. A 2nd-order regressor set U generated by this scheme is not only 
much smaller than an exhaustive pair-set, but also biologically more plausible. Note that it 
is possible to include other sets of SNP pairs from other resources in our candidate set. For 
example, in our experiments, we also added SNP pairs that passed two-locus epistasis test 
with p-value < 10 -5 into the set U. 

After finding the candidate SNP pairs, we generate the group of SNPs or interacting SNP 
pairs in two steps. In the first step, we find highly interconnected subgraphs (or clusters) from 
the genetic interaction network using any graph clustering algorithms. In our experiments, 
we used MCODE algorithm 2j for clustering the network. In the second step, we group all 
the SNPs or SNP pairs that are linked to the genes in a cluster. We linked the genes and 
SNPs based on physical locations in the genome. For example, if a SNP is located nearby 
a gene within a certain distance (e.g. <500bp), they are linked together. Finally, we define 
individual SNPs in the mth group as g m e Q and SNP pairs in the mth group as \ m e L. 
We then look for associations between inputs/input-pairs and outputs via Eq. ( jl~5f) : 

K N / J \ 2 

2 E E U - E &i - E ( 15a ) 

fc=i i=i y j=l (r,s)eU / 



mm 
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K J 



AiEEi^i ( 15b ) 
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k / \g\ , \c\ 



+ x * E E x / E w + EJE ^ 2 1 ( 15c ) 

k=l \ m=l V j'Ggm m=l V (r,s)el m 



j \H\ , \H\ 



\j'=l m=l y fcGh m (r,s)eUm=l y fceh m 7 

if 

fc=l (r,s)eU 

where ^ is the set of input groups for marginal terms and L is the set of input groups for 
pairwise interaction terms. Here, we use two tuning parameters for L\ penalty depending on 
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whether a covariate is modeling an individual effect (Ai) or interaction effect (A4) because 
they might need different levels of sparsity. Note that this problem is identical to Eq. 
§5§ if we treat interaction terms additional covariates, and hence our optimization 

method presented in Section H] is applicable to Eq. ( TT5l) . However, Eq. ( TT51) will be more 
computationally expensive than Eq. ([5]) since Eq. ( Tl5l) has a larger number of covariates in 
B including both marginal and interaction terms and additional tuning parameter A4. 

6. SIMULATION STUDY 
In this section we validate our proposed method using simulated genome/phenome datasets, 
and examine the effects of simultaneous use of input and output structures on the detection 
of true non-zero regression coefficients. We also evaluate the speed and the performance 
of our optimization method for support recovery in comparison to two other alternative 
methods. For the comparison of optimization methods, we selected smoothing proximal 



gradient method 



19J since both methods are in principle able 



and the union of supports 
to use input/output structures and handle overlapping groups. 

The simulated datasets with J = 120, K = 80, and N = 100 are generated as follows. For 
generating X, we first selected 60 input covariates from a uniform distribution over {0, 1} 
which indicates major or minor genotype. We then simulated 60 pairwise interaction terms 
(x l j x Xj,) by randomly selecting input-pairs from the 60 covariates mentioned above. Pooling 
the 60 marginal terms and 60 pairwise interaction terms resulted in a input space of 120 
dimensions. We also defined input and output groups as follows (for the sake of illustration 
and comprehension convenience, here our input and output groups correspond to variables 
to be jointly selected rather than jointly shrunk, the shrinkage penalty in our regression loss 
can be defined on the complements of these groups): 



5, . . . , 9, 10, . . . , 15, . . . , 25, . . . , 29, 30, 31, 32, . . . , 37, . . . , 50, . . . , 54, 55, 56, 57, . . . , 60, . . . , 75, . . . , 80, . . . , 87, . . . , 94, . . . , 104, .... 109, 110, 111, ... , 116 



en 

»5 



1, . . . , 4, 5, . . . , 10, . . . , 12, . . . , 17, 18, 19, 20, . . . , 25, . . . , 46, . . . , 56, . . . , 63, . . . , 70, . . . , 75, . . . , 80, 



h2 h4 he h7 
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where the numbers within a bracket represent the indices of inputs or outputs for an input 
group g t , t = 1, . . . , 10, or an output group h Q , o = 1, . . . , 7. The inputs and outputs 
which did not belong to any groups were in a group by itself. We then simulated B, i.e, the 
ground truth that we want to discover. We selected non-zero coefficients so that B includes 
various cases, e.g., overlap between input and output groups, overlap within input groups, 
and overlap within output groups. Figure |2Ja) shows the simulated B £ jj80xi2o w h ere 
non-zero coefficients are represented by black blocks. Given X and B, we generated K = 80 
outputs by Y = BX + E, E ~ A/"(0, 1). We generated 20 datasets and optimized Eq. (JSJ) 
using the three methods. We report the average performance using precision recall curves. 




(a) (b) (c) (d) 

Figure 2: An example of simulation results with \P J k \ — 2, N — 100, J = 120, and K = 80. 
(a) True regression coefficient matrix. Estimated B by SIOL (b) with both input and output 
structures (c) with only input structure, and (d) with only output structure. In (b-d), we 
show the normalized values of | B | . 

6.1 Evaluation of the Effects of Using Input and Output Structures 

We first investigate the effects of using both input and output structures on the performance 
of our model. Here we applied our optimization method (HiGT) to the following three 
models with different use of structural information: 

1. Use of both input and output structures (Eq. (|5a|) + Eq. ()5b[) + Eq. (|5cl) ) 

2. Use of input structures (Eq. (|5aj) + Eq. (|5b])) 

3. Use of output structures (Eq. (|5a|) + Eq. ([3c])) 
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We then observed how the use of input / output structures affect the recovery of the true 
non-zero coefficients and the prediction error. In Figure [2], we visualize the examples of 
estimated B by the three different models. Figure [2](b) shows that the model with input and 
output structure successfully recovered true regression coefficients in Figure [2(a). However, 
as shown in Figure [2]^ c-d), the models with either input or output structure were less effective 
to suppress noisy signals, which resulted in many false positives. 




Recall Recall Recall 

(a) (b) (c) 

Figure 3: Precision recall curves on the recovery of true non-zero coefficients due to SIOL with 
both input and output structures (input/output struct), regression with only input structure 
(input struct), and with only output structure (output struct), under three different signal 
strengths of true regression coefficients, (a) (3 3 k = 0.4, (b) (3 J k = 1, and (c) f3 3 k = 2. The 
simulated data were generated with N = 100, J = 120, and K = 80. 

Figure [3] shows the precision recall curves on the recovery of true non-zero coefficients by 
changing the threshold r for choosing relevant covariates > r), under different signal 
strengths of 0.4, 1 and 2. For all signal strengths, the model with input/output structures 
significantly outperformed the other models with either input or output structure. The most 
interesting result is that when the signal strength was very small such as 0.4, our model still 
achieved good performance by taking advantage of both structural information. 

We also compare the prediction errors on our validation data with 280 (20 x 14) sam- 
ples (each dataset had 14 samples for validation). For computing the prediction error, we 
first selected non-zero coefficients, and then recomputed the coefficients of those selected 
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(a) (b) (c) 

Figure 4: Comparison of the prediction error of SIOL (input/output struct), with regression 

under only input structure (input struct), on only output structure (output struct), (a) 

PI = 0.4, (b) PI = 1, (c) Pi = 2. 

covariates using linear regression without shrinkage penalty. Using the unbiased coefficients 
of the chosen covariates, we measured the prediction error for our validation data. Figure 
|4] shows the prediction error under different signal strengths ranging from 0.4 to 2. For all 
signal strengths, we obtained significantly better prediction error using both input and out- 
put structures. When the signal strength was large such as 1 or 2, the use of both input and 
output structures was especially beneficial for reducing the prediction error since it helped 
the model to find most of the true covariates relevant to the outputs. 

The effects of the size of input and output groups Figure (E^-EH) demonstrates the 
results on simulated datasets with different size of input and output groups. For all group 
sizes, our method significantly improved the performance by effectively taking advantage of 
both input and output groups. 

6.2 Comparison of HiGT to Alternative Optimization Methods 

In this section, we compare the accuracy and speed of our optimization method (HiGT) 
with those of the two alternative methods including smoothing proximal gradient method 



(SPG) [8[] and union of supports 



3- 



Both alternatives can handle overlapping groups. 
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Figure 5: Precision recall curves on the recovery of true non-zero coefficients for SIOL 
(input/output struct), regression with only input structure, and with only output structure, 
under two different sizes of input and output groups, (a) \g\ G {2,3}, (b) \g\ = 5, Wg G Q, 
and (c) \h\ = 5, and (d) \h\ = 40, Wh G %. We fixed the size of output groups for (a,b) 
(\h\ = 10), and fixed the size of input groups for (c,d) (\g\ = 5). The simulated data were 
generated with fi{ = 0.5, N = 100, and J = 120. 

Specifically, the smoothing proximal gradient method is developed to efficiently deal with 
overlapping group lasso penalty and graph-guided fusion penalty using an approximation 
approach. However, it may be inappropriate for our model since the maximum gap between 
the approximated penalty and the exact penalty is proportional to the total number of 
groups R, where R = J\H\ + i£"|£7|- Thus, when dealing with high dimensional data (e.g 
J ~ 500, 000) such as genome data, the gap will be large, and the approximation method 
can be severely affected. On the other hand, "union of supports" finds the support of B from 
the union support of overlapping groups. To obtain the union of supports, input variables 
are duplicated to convert the penalty with overlap into the one with disjoint groups, and a 
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standard optimization technique for group lasso 52] can be applied. One of disadvantages 
of union of supports is that the number of duplicated input variables increases dramatically 
when we have a large number of overlapping groups. In our experiment, we considered all 
possible combinations of overlapping input and output groups, and used a coordinate descent 
algorithm for sparse group lasso 
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Figure 6: Precision recall curves on the recovery of true non-zero coefficients using the SIOL 
model via HiGT, smoothing proximal gradient method (SPG), and union of supports for 
optimization. Three different model sizes determined by the number of input variables were 
tested (due to high computational cost, results of union-of-support are only available for the 
smallest problem sizes tested): (a) J = 30, (b) J = 400, and (c) J = 600. The simulated 
data were generated with {3 3 k = 2, N = 100, and K = 20. 

Figure [6] shows the precision recall curves on the recovery of true non-zero coefficients 
under the SIOL model using the three optimization methods. The size of the problem is 
controlled by increasing number of input variables (from 30 to 600). The simulated data set 
used here was identical to the data in Section IBTTl except that we used 20 outputs (y6i, • • • , yso) 
and different number of input variables. One can see that our method outperforms the other 
alternatives for all configurations. Our method and smoothing proximal gradient method 
showed similar performance when the input variable is small (J = 30) but as J increases, 
our method significantly performed better than SPG. It is consistent with our claim for the 
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maximum gap between the approximated penalty and the exact penalty which is related to 
the number of groups. Union of supports did not work well even when the number of input 
variables is small (J = 30) since the actual number of input variables considered was very 
large due to the duplicated covariates, which severely degraded the performance. 
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(a) (b) 
Figure 7: Time complexity of HiGT, SPG, and union of supports. All three methods used 

both input and output groups, (a) Computational time with different number of samples, (b) 

computational time with different number of inputs. We used the same tuning parameters 

for all the methods (Ai = 0.01, A2 = A3 = 0.1). We did not report the times for the small 

number of samples and inputs for our method and SPG since I/O latency was dominant. 

We also compared the speed of our method with the two alternatives of union of supports 
with all possible combinations of input and output groups and SPG that considered both 
input and output groups. Figure[7](a,b) show that our method converged faster than the other 
competitors, and was significantly more scalable than the two alternatives. Union of supports 
was very slow compared to our method and SPG because of the large number of duplicated 
input variables. Our experimental results confirms that our optimization technique is not 
only accurate but also fast, which can be explained by the use of DAG and the simple forms 
of optimality checks. 
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7. ANALYSIS OF YEAST EQTL DATASET 
We apply our method to the budding yeast (Saccharomyces cerevisiae) data [4j with 1,260 
unique SNPs (out of 2,956 SNPs) and the observed gene-expression levels of 5,637 genes. As 



101 ] with stringent 



network prior knowledge, we used genetic interaction network reported in 
cutoff to construct the set of candidates of SNP pairs U. We follow the procedure in section 
[5] to make U with an additional set of significant SNP pairs with p- value < 10 -5 computed 
from two-locus epistasis test. When determining the set U, we assumed that a SNP is linked 

e choice 



43|. As 



to a gene if the distance between them is less than 500bp. We consider it a reasonab 
for cis-effect as the size of intergene regions for S. cerevisiae is 515bp on average 
a result, we included 982 interaction terms from the interaction network in X with 1,260 
individual SNPs. The number of SNP pairs from two-locus epistasis test varied depending 
on the trait. For generating input structures, we processed the network data as follows. We 
started with genetic interaction data which include 74,984 interactions between gene pairs. 
We then extracted genetic interactions with low p-values (<0.001). Given 44,056 significant 
interactions, using MCODE clustering algorithm, we found 55 gene clusters. Using the gene 
clusters, we generated the groups of individual SNPs and pairs of SNPs according to the 
scheme in section |5j For generating output structures, we applied hierarchical clustering to 
the yeast gene expression data with cutoff 0.8, resulting in 2,233 trait clusters. 

Marginal Effects in Yeast eQTL dataset We briefly demonstrate the effects of in- 
put/output structures on the detection of eQTLs with marginal effects. In general, the 
association results for marginal effects by our method, lasso and single SNP analysis (the 
later two are standard methods in contemporary GWA mapping that use no structural infor- 
mation, and hence included for comparison) showed similar patterns for strong associations. 
However, we observed differences for SNPs with small or medium sized signals. For exam- 
ple, our results had fewer nonzero regression coefficients compared to lasso. One possible 
explanation would be that the grouping effects induced by our model with input/output 
structures might have removed false predictions with small or medium sized effects. To 
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Figure 8: Manhattan plot for association between (YER160C and YJR029W) and SNPs on 
chromosome 7. The two genes YER160C and YJR029W share the same GO category "trans- 
position" . Our method detected SNPs which affect both two genes in this region. However, 
single SNP analysis did not find any associated SNPs and lasso found SNPs associated only 
with YER160C in this region. Graph were generated using the GenAMap software 
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illustrate eQTLs with marginal effects, we show some examples of association SNPs using 



GenAMap . Figure |8] demonstrates a Manhattan plot on chromosome 7 for two genes 
including YER160C and YJR029W. Both genes have the same GO category "transposition" . 
As both genes share the same GO category, it is likely that they are affected by the same 
SNPs if there exist any association SNPs for both genes. In our results, we could see that 
the same SNPs on chromosome 7 are associated with both genes as shown in Figure El How- 
ever, single SNP analysis did not find any significant association SNPs in the region. Lasso 
detected association SNPs in the region but they were associated with only YER160C rather 
than both of them (lasso plot is not shown to avoid cluttered figure). This observation is in- 
teresting since it supports that our method can effectively detect the SNPs jointly associated 
with the gene traits by taking advantage of structural information. 

Epistatic Effects in Yeast eQTL dataset Now we show the benefits of using the in- 
put/output structures for detecting interaction effects among SNPs by comparing the results 



of our method to those of two-locus epistasis test performed by PLINK 37J which uses no 
structural information. Specifically, we compare the hotspots with interaction effects (i.e. 
SNP pairs that affect a large number of gene traits) which are identified by both methods. 
Recall that two-locus epistasis test is the most widely used statistical technique for detect- 
ing interaction effects in genome-wide association studies, which computes the significance 
of interactions by comparing between the null model with only two marginal effects and the 
alternative model with two marginal effects and their interaction effect. In the following 
analysis, we discarded all SNP pairs if the correlation coefficient between the pairs > 0.5 to 
avoid trivial interaction effects. 

We first identified the most significant hotspots that affect more than 100 gene traits. To 
make sure that we include only significant interactions, we considered interaction terms if 
their absolute value of regression coefficients are > 0.05. For the results of two-locus epistasis 
test, we considered all SNP pairs with p-value < 10 -5 . Figure [9]^a,b) show the hotspots 
found by our method and two- locus epistasis test. The rings in the figure represent the yeast 
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Figure 9: Hotspots with interaction effects identified by (a) our method and (b) two-locus 
epistasis test. This figure represents the yeast genome in a circular format. In clockwise 
direction, from the top of the circles, we show 16 chromosomes, which are separated with 
space and different colors. Lines indicate interaction effects between two connected locations 
in the genome. Thickness of the lines is proportional to the number of traits affected by 
the interaction effects. Here we show interaction effects which influence more than 100 gene 
traits. The hotspots for (a) are represented in Table [H In (b), two SNP pairs are found 
including chrl6:718892-chrl6:890898 (affected genes are enriched with the GO category of 
ribosome biogenesis with corrected p- value 1.6 x 10 -36 ), and chr8:56246-chr9:362631 (affected 
genes are enriched with the GO category of vacuolar protein catabolic process with corrected 



p- value 1.6 x 10 ). This figure was generated using Circos software 



25). 
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Table 1: Hotspots of SNP pairs having epistatic effects in yeast identified by our method. 



Hot spot 


SNPf 


SNP2 


Number of 


GO category of 


Corrected p- value of 


label 


location 


location 


affected traits 


affected traits 


GO category 


f 


chrl:l 54328 


chr5:350744 


455 


ribosome biogenesis 


1.2 x 10~ 36 


2 


chrf0:380085 


chrf5:f70945 


195 


ribosome biogenesis 


1.6 x 10~ 12 


3 


chrf0:380085 


chrf5:f75594 


185 


ribosome biogenesis 


4.1 x 10~ 12 


4 


chr5:222998 


chrf5:f08577 


170 


response to temperature stimulus 


2.9 x 10~ 6 


5 


chrff:388373 


chrf3:64970 


155 


regulation of translation 


1.8 x 10~ 32 


6 


chr2:4990f2 


chrf5:5f9764 


145 


vacuolar protein catabolic process 


1.4 x 10~ 7 


7 


chrf:4f483 


chr3:643ff 


130 






8 


chr7:f4f949 


chr9:277908 


125 






9 


chr3:643ff 


chr7:3f2740 


115 


glycoprotein metabolic process 


1.5 x 10" 4 


fO 


chrf2:957f08 


chrf5:f70945 


110 


vacuolar protein catabolic process 


7.8 x 10~ 16 


ff 


chr4:864542 


chrf3:64970 


105 


ribonucleoprotein complex biogenesis 


3.7 x 10" 6 



genome from chromosome 1 (located at the top of each circle) to 16 clockwise, and the lines 
show interactions between the two genomic locations at both ends. One can see that our 
method detected 11 hotspots but two- locus epistasis test found only two significant hotspots 
with interaction effects. This observation shows that our method can find more significant 
hotspots with improved statistical power due to the use of the input /output structures. In 
Table [U we summarized the hotspots found by our method. It turns out that our findings 
are also biologically interesting (e.g. 9 out of 11 hotspots showed GO enrichment). Notably, 
hotspot 1 (epistatic interaction between chrl:154328 and chr5:350744) affects 455 genes which 
are enriched with the GO category of ribosome biogenesis with a significant corrected p- value 



< 10 (multiple testing correction is performed by false discovery rate 



29|). This SNP pair 



was included in our candidates from the genetic interaction network. There is a significant 



10|, and both genes 



genetic interaction between NUP60 and RAD51 with p-value 3 x 10~ 7 
are located at chrl:152257-153877 and chr5:349975-351178, respectively. As both SNPs are 
closely located to NUP60 and RAD51 (within 500bp), it is reasonable to hypothesize that two 
SNPs at chrl: 154328 and chr5:350744 affected the two genes, and their genetic interaction 
in turn acted on a large number of genes related to ribosome biogenesis. 

To provide additional biological insights, we further investigated the mechanism of this 
significant SNP-SNP interaction. From literature survey, RAD51 (RADiation sensitive) is 



a strand exchange protein involved in DNA repair system 42], and NUP60 (NUclear Pore) 



121. Also, it has 



is the subunit of unclear pore complex involved in nuclear export system 
been reported that yeast cells are excessively sensitive to DNA damaging agents if there exist 
mutations in NUP60 



33] . In our results, we also found out that the SNP close to NUP60 did 
not have significant marginal effects, and the SNP in RAD51 had marginal effects. According 
to these facts, it would be possible to hypothesize the following. When there is no mutation 
in RAD51, the point mutation in NUP60 cannot affect other traits since the single mutation 
is not strong enough and if there exist DNA damaging agents in the environment, DNA 
repair system would be able to handle them. However, when there exist point mutations in 
RAD51, DNA damaging agents would severely harm yeast cells with the point mutation in 
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NUP60 since DNA repair system might not work properly due to the mutation in RAD51 
(recall that the SNP in RAD51 had marginal effects). As a result, both mutations in NUP60 
and RAD51 could make a large impact on many gene traits. 

chrl chrl 




chrlo cW 9 

(a) 



chrlo cW 9 



(b) 



Figure 10: Hotspots with interaction effects identified by (a) our method and (b) two-locus 



epistasis test by PLINK. Here we show epistatic interactions w 
gene traits. This figure was generated using Circos software 



rich influence more than 10 



25]. 



Furthermore, we looked at the hotspots which affect > 10 gene traits. Figure [TOT a.b) 
show epistatic interactions identified by our method and two- locus epistasis test, respectively. 
In this figure, we show significant interactions with regression coefficient cutoff > 0.1 for our 
method, and p-value cutoff < 10~ 6 for two-locus epistasis test. These cutoffs are arbitrarily 
chosen to make the number of hotspots found by both methods similar. Surprisingly, two 
methods showed very different hotspots with epistatic interactions. Figure [TOT a) was very 
similar to Figure ED^a) but in Figure [TDT b). several hotspots emerged which were absent 
in Figure Mjo). We will analyze these hotspots in two ways. First we will look at the 
hotspots with epistatic effects which appeared in both Figure [TOT a) and (b). Then we will 
investigate the differences between the two results. First, we observed that both methods 
found significant epistatic effects between chromosome 1 and 5. Recall that in our previous 



36 



1 

0.8 

CM 



i 0-4 

s— 
s— 

o 

o 

0.2 


0.2 0.4 0.6 0.8 1 
Correlation for SNP1 

Figure 11: The scatter plot for illustrating the correlation between our epistatic hotspot 
1 (chrl:154328-chr5:350744) and significant SNP pairs close to hotspot 1 detected by two- 
locus epistasis test (p-value < 10 -6 and the distance between the pairs of SNPs and hotspot 
1 is within < 50kb). Each dot represents a SNP pair (SNP1, SNP2) found by two-locus 
epistasis test, and x-axis represents the correlation between SNP1 and chr 1:154328 and y- 
axis represents the correlation between SNP2 and chr5:350744. Each dot was perturbed by 
a small amount of random noise to avoid overlapping of the dots. 
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analysis of the hotspots, this interaction was discussed (see hotspot 1 in Tabled]). Among 
all significant SNP pairs found by two- locus epistasis test, there was no identical SNP pair 
to hotspot 1 but there were 30 SNP pairs close to it (within < 50kb). Also, it turns out 
that these 30 SNP pairs had very strong correlation with hotspot 1. In Figure [HJ we show 
scatter plot to illustrate the strong correlations between hotspot 1 and these 30 SNP pairs. 
More interestingly, the total number of affected traits by these 30 SNP pairs was 416, and it 
is very similar to 455 that is the number of affected genes by hotspot 1. According to these 
facts and our previous analysis for the mechanism of hotspot 1, it seems that hotspot 1 is 
truly significant, and two-locus epistasis test found significant SNP pairs that are close to 
the true location but failed to locate the exact location of hotspot 1. It supports that our 
algorithm could find such a significant hotspot affecting > 400 genes by detecting exact SNP 
pairs. However, two-locus epistasis test was unable to locate many hotspots affecting a large 
number of traits due to insufficient statistical power. Second, we investigated the differences 
between the two results in Figure fTOTa.b). As we cannot report all the results in the paper, 
we focused on a SNP pair (chrl0:87113-chrl5:141621) in Figure fTDT a). and another SNP pair 
(chr8:63314-chr9:362631) in Figure ITOT b). Figure fT2T a.b) show the average gene expression 
levels for each SNP pair. In this figure, x-axis represents the genotype G {0, 1} which is 
the multiplication of two SNPs (SNP1 x SNP2, where SNP1, SNP2 e {0,1}), and y-axis 
represents the average gene expression levels of individuals with given genotype. Each line 
in Figure IT27 a.b) shows how the average gene expression level varies as the genotype changes 
from to 1 for each trait affected by the SNP pairs with error bars of one standard deviation. 
Interestingly, in Figure IT27 a). we could see that there is a consistent pattern, where for most 
gene traits, the expression levels decreased as the genotype changed from to 1. However, as 
shown in Figure [T27 b). for the SNP pair found by two-locus epistasis test, we could not find 
such a coherent pattern. It demonstrates the differences between our method and two-locus 
epistasis test. As our model borrows information across input and output group structures, 
we could find consistent gene expression patterns for the SNP pair. On the other hand, two- 
locus pairwise test analyzed each SNP pair separately, and each trait affected by the SNP 
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pair showed different gene expression patterns with different standard deviations. Thus, it 
seems that our method can provide interesting biological insights in terms of gene expression 
patterns in addition to the statistical significance. 



0.6 r 2 
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Figure 12: Variations of gene expression levels w.r.t. to the genotypes of (a) a SNP pair 

(chrl0:87113-chrl5:141621) found by our method, and (b) a different SNP pair (chr8:63314- 

chr9:362631) found by two-locus epistasis test. Here, x-axis represents genotype doses and 

y-axis shows the average expression levels of the multiple genes (denoted by multiple vertical 

lines) affected by the corresponding SNP pairs. A small noise was added to the genotypes 

as offsets to avoid overlapping of the error bars. 

8. DISCUSSIONS 

In this paper, we presented jointly structured input-output lasso to simultaneously take ad- 
vantage of both input and output structures. We also presented an efficient optimization 
technique for solving our multi-task regression model with structured sparsity. Our experi- 
ments confirmed that our model is able to significantly improve the accuracy for detecting 
true non-zero coefficients using both input and output structures. Furthermore, we demon- 
strated that our optimization method is faster and more accurate than the other competitors. 
In our analysis of yeast eQTL datasets, we identified important pairs of eQTL hotspots that 
potentially interact with each other. 
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Prior knowledge about input and output structures In practice, it is important 
to generate reliable input and output groups to maximize the benefits of the structural 
information. In our experiments, we showed that yeast genetic interaction networks can be 
used as prior knowledge to define input and output structures. However, such reliable prior 
knowledge cannot be easily attained when we deal with human eQTL datasets. Instead, we 
have a variety of resources for human genomes including protein-protein interaction networks 
and pathway database. Generating reliable input and output structures exploiting multiple 
resources would be essential for the successful discovery of human eQTLs. 

Comparison between HiGT and other optimization algorithms Recently, an ac- 



tive set algorithm 



20( has been proposed developed for variable selection with structured 



sparsity, which can potentially be used for estimating the SIOL model. We observe two 
key differences between our method and the active set algorithm 2(|. First, the active 
set algorithm incrementally increases active sets by searching available non-zero patterns, 
hence one can see that it is a "bottom-up" approach. On the other hand, our method 
adopts "top-down" approach where irrelevant covariates are discarded as we walk through 
the DAG. Second, our algorithm guarantees to search all zero patterns while the active set 
algorithm needs a heuristic to select candidate non-zero patterns. When B is sparse, our 
algorithm is still very fast by taking advantage of the structures of DAG. However, when 
B is not sparse, our algorithm needs to search a large number of zero patterns and update 
many non-zero coefficients but the active set algorithm still does not need to update many 
non-zero coefficients. Hence, in such a non-sparse case, the active set algorithm may have 
an advantage over our optimization method. Other alternative optimization methods for 



SIOL include MM (majorize/minimize) algorithm 26] and generalized stage-wise lasso 55]. 
However, these methods did not work well for SIOL as the approximated penalty by MM al- 
gorithm and the greedy procedure of generalized stage- wise lasso were incapable of efficiently 
inducing complex sparsity patterns. 
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Future work One promising research direction would be to systematically estimate the 
significance of the covariates that we found. For example, computing p- values of our results 
would be helpful to control the false discovery rate. To handle both sparse and non-sparse 
B, it would be also interesting to develop an optimization method for our model that can 
take advantage of both "bottom-up" and "top-down" strategies. For example, we can select 
variables using "bottom- up" approach and discard irrelevant variables using "top-down" 
approach alternatively in a single framework. Finally, we are interested in applying our 
method to human disease datasets. In that case, the extension of our work to handle case- 
control studies and finding reliable structural information will be necessary. 
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